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Abstract 

We present a Lohner type algorithm for the computation of rigorous 
bounds for Solutions of ordinary differential equations and its derivatives 
with respect to initial conditions up to arbitrary order. As an applica- 
tion we prove the existence of multiple invariant fori around some elliptic 
periodic orbits for the pendulum equation with periodic forcing and for 
Michelson System. 



1 Introduction 

This paper is a sequel to |ZJ. We present here a Lohner- type algorithm for 
computation of rigorous enclosures of partial derivatives with respect to initial 
conditions up to an arbitrary order r of the flow induced by an autonomous 
ODE, hence the name C-Lohner algorithm. Let r be a positive integer, then 
by C-algorithm we will mean the routine which gives rigorous estimates for 
partial derivatives with respect to initial conditions up to an order r and C- 
computations we mean an application of an C'-algorithm. 

Our main motivation for the developmcnt of C-algorithm was a desire to 
provide a tool, which will considerably extend the possibilities of Computer 
assisted proofs in the dynamics of ODEs. Till now most of such proofs have used 
topological conditions (see for example [HZHTl IMM[ IGZ[ IZT] ) and additionally 
conditions on the first derivatives with respect to initial conditions (see for 
example [RNSl [Tl IWiTl IWZl [KZ] ) . hence it required C- and C^-computations, 
respectively. The spectrum of problems treated includes the questions of the 
existence of periodic orbits and their local uniqueness, the existence of symbolic 
dynamics, the existence of hyperbolic invariants sets, the existence of homo- 
and heteroclinic orbits. To treat other phenomena, likc bifurcations of periodic 
orbits, the route to chaos, invariant tori through KAM theory one needs the 
knowledge of partial derivatives with respect to initial conditions of higher order. 

In principle, one can think that a good rigorous ODE solvcr should be en- 
ough. Namely, to compute the partial derivatives of the flow induced by 

x' = f{x), X £ R" (1) 
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it is enough to rigorously integrate a System of variational equations obtained by 
a formal difFerentiation of ([T]) with respect to the initial conditions. For example 
for r = 2 we have the foUowing System 

= /(i). (2) 

s— 1 

s.r—1 s—1 

with the initial conditions 

x{0)^xo, V{0)=ld, H,jk{0)^0, i,j,k^l,...,n. 

It is well known that if by Lp{t,XQ) we denote the (local) flow induced by ([T]), 
then 



dip. 
dxj 

dxjdxk 



Analogous Statements are true for higher order partial derivatives with respect 
to initial conditions. 

It turns out that a straightforward application of a rigorous ODE solver to 
the System of variational Equations (ßHll) is very ineflicient. Namely, it totally 
ignores the structure of the System and leads to a very poor Performance and 
unnecessary long computation times (see Section |4J|) . 

Our algorithm is a modification of the Lohner algorithm [Lo| . which takes 
into account the structure of variational Equations ([2H11) . Basically it consists of 
the Taylor method, a heuristic routine for a priori bounds for Solution of ([2H4]) 
during a time step and a Lohner-type control of the wrapping effect, which is 
done separately for x and partial derivatives with respect initial conditions (the 
variables V and H in ()3|4p ). The Taylor method is realized using the automa- 
tic differentiation fRaj and the algorithms for computation of compositions of 
multivariate Taylor series. 

The proposed algorithm has been successfuUy applied in [HNW| to the Mi- 
chelson System |Mi] , where a Computer assisted proof of the existence of a cocoon 
bifurcation was presented. Some parts of this proof required C^-computations. 

In the present paper in Section[8]we show an application of our algorithm to 
pendulum equation with periodic forcing and the Michelson System. We used 
it to compute rigorous bounds for the coefhcients of some normal forms up to 
Order five, which enabled us to prove the existence of invariant tori around some 
elliptic periodic orbits in these Systems using KAM theorem for twist maps on 
the plane. These proofs required and computations. 
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2 Basic definitions 



To effectively deal with thc formulas involving partial derivatives we will use ex- 
tensively a notation of multiindices, multipointers and submultipointers throu- 
ghout the paper. 

As an motivation let us consider the formula for tlie partial derivatives of 
the composition of maps. Assume 5 : R" ^ K" and / : ^ M are of class C^. 

We have 

d^{f°9) ^ d'^^f dgk dgr dg^ ^ df d^gk 

dxidxjdxc ^ dxkdxrdxs dxi dxn dxc ^ dxk dxidx^dxc 

k,r,s—l k—1 

^ d'^f f d^gk dgr _^ dgk d'^gr _^ d'^gk dgr 

dxkdxr \dxidxc dxj dxi dxjdxc dxidxj dxc 

To thc Operator ^ we can in a unique way assign a multipointer, 

which is a nondecreasing sequence of integers (ji,j2,j3); such that {11,12, «3} = 
{ji) 72,^3} • A submultipointer is a multipointer, which is a part of a longer mul- 
tipointer, for examplc c)(i 3) = (i,c). One observes, that submultipointers 
appear at several places in the above formula. 

A multiindex is an dement of a e N". It is another way to represent various 
partial derivatives. Thc coefficient teils us how many times to differentiate 
a function with respect to the i-th variable. Obviously, we have one-to-one 
correspondence between multipointers and multiindices. 

2.1 Multiindices 

By N we will denote the set of nonnegative integers, i.e. N = {0, 1, 2, . . .}. 

Definition 1 An element r € N" will he called a multiindex. 

For a sequence a = (ai, . . . , «„) e N" and a vector x = {xi,. . . , a;„) e M" we 
set 

1. |q;| = ai + • • • + a„ 

2. a! = ai! • «2! • • • 

3. x°' = {xi^ , . . . , a;^" ) 
By e" e N" we will denote 



er = (0,0,...,0, 1 ,0,...,0,0). 

We will drop the index n (the dimension) in the Symbol e" when it is obvious 
from the context. 

Put N^' := {a e N" : \a\ =p}. 

For S = {61, . . . , e X • • • X we set 
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1- l^l = E-=il^il 

Let / = (/i, . . . , fm) : M" ^ be sufficiently smooth. For a e N" we set 

1. D" fi = — 

2. D^f={D"fi,D"f2,...,D"fm) 

For a function / : M x K" ^ M" by D°'fi{t,x) we will denote D'^fi{t, ■){x) and 
similarly 

x) = {D'^hit, x),..., D'^Uit, x)). 
This Convention means that always acts on a;- variables. 

2.2 Multipointers 

For a fixed n > and p > we define 

:= {(ai, 02, . . . ,ap) & W : l < ai < ■ ■ ■ < ap < n) 

OO 

M = N'' := 

Definition 2 An dement o/A^" will he called a multipointer. 
Remark 3 A function 

A:A/';9(ai,...,ap)^f^e;\eN;' (5) 

i=l 

is a bijection. 

Let / = (/i, . . . , fm) : M" ^ be a sufficiently smooth. For a€J\f^ we set 

OXai • • • OXap 
2. DJ:= {Dafl,...,Dafm) 

For a function / ; R x M" ^ K" by Dafi(t,x) we will denote Dafi{t, ■){x). In 
the light of thc above notations Daf = D^^"^ f. 

For a= (oi, 02, . . . , o„) G and & = (&i, 62, • ■ • , ^n) € we define 

o + 6 = (ßi + 61, . . . , o„ + 6„) G N;+^. 

For a G AAp and /? G Af^ wc define 

a + /3 = A-i (A(a) + A(/3)) G A/;%. 
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By < WC will denote a linear order (lexicographical order) in A/" defined in 
the foUowing way. For a e Afp and b € Afg 

( < b) -. i cithcr 3i,i < p,i < q, ai < bi and üj = bj for j <i 
~ \orp<q and Oj = öj for i = 1, . . . , p. 

Definition 4 For k < p we set 

M^ik) := ...,ök)e {Mn'' ■.Si<---<Sk,ö^ + -- ■ + Sk = (1, 2, . . . ,p)} (7) 

We will use Af^{k) extensively in the next section. Its will be used to label 
terms in D"fi{(fi{t,x)). Observe that for p > 

Af^il) = {{l,2,...,p)} 
Äf'>{p) = {{{l),{2),...,{p))} 

One can construct all elements of J\f^{k) using the following recursive procedure. 
From the definition of J\f^{k) it foUows that if . Sm-i) € J^P~'^{m. - 1) 
then {öl, . . . ,6m-i, {p)) G N^{m) (notice that order is preserved). Similarly, if 
('5i,...,<5™) e A/'P-i(m) then 

{5i, . . . .,6s-u5s + {p),5s+i, . ■ . ,5^) &J\fP{m) 

and again order of elements is preserved. Hence, for p > 2 and 1 < < p we 
have AfP{k) = AyjB where 

A = {{Su ök-i, (p)) : . . . € f^^-'ik ~ 1)} 

B=\J {{Su...,6,-uS, + {p),Ss+i,...,6k):{Su...,dk)€J^P-\k)} 

s=l 

and the sets A and B are disjoint. 

Another way to generate all elements ofJ\fP{k) can be described as foUows 

• decompose the set {1,2, into k nonempty and disjoints sets Aj, 
i = 1, . . . ,k 

• we sort each Aj and permute Aj's to obtain min(Ai) < min(A2) < • • • < 

min(Afe) 

• we dcfinc öi to be an ordered set consisting of all elements of Aj for 
i = 1, . . . ,k 

Definition 5 For an arbitrary a € Afp and ö € Af^ such that k < p we define a 
submultipointer as G Af^ by {as)i = as^ for i = 1, . . . ,k, which can be expressed 
using A as follows 
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3 Equations for variations 



Consider an ODE x' = /(x) where / is C^+^ Let : M x M"^M" be a local 
dynamical System induced by x' = /(x). It is well known, that Lp e and one 
can derive the equations for partial derivatives of </? by differentiating equation 
^(tjx) = f{ip{t,x)) with respect to the initial condition x. As a result we 
obtain a System of so-called equations for variations, whose size depends on the 
Order r of partial derivatives we intcnd to compute. An example of such System 
for r = 2 is given by ([IHl]) with initial conditions given by ([5]). 

The goal of this section is to write the equations for variations in a compact 
form using multipointers and multiindices, which allows us to take into account 
the symmetries of partial derivatives, 

Lemma 6 Assume f G C^^^ and let (p : M. x ]R"-<^M" be a local dynamical 
System induced by x' — f{x). Then for a G J\fp such that p < r holds 

j p 71 k 

k=l ii,...,ifc = l (<5i,...,Äfc)e7VP(fc) 3 = 1 

for i = 1, . . . ,n. 

Proof: In the proof the functions D'^^^ ^^'^ /j are always evaluated at 
(p{t, x), and various partial derivatives of ip are always evaluated at (f, x), there- 
fore the arguments will be always dropped to simplify formulae. We prove the 
lemma by induction on p — \a\. If p = 1 then a = (c) for some c G {1, . . . , n} 
and ^ becomes 

s— 1 s— 1 

Assume (|9]) holds true ior p — 1, p > 1. Let us fix a G Af^. We have 
a — b + (c), where b = (oi, . . . , Op^i) G Äfp^i and c = Op. Since ([5]) is satisfied 
for p — 1, therefore we have 

d f d 

p— 1 ri fe ^ 

E E D^f^ E n^s^. 

, fc=l ii,...,ifc = l (<5i,...,5fc)eA^p-i(fc) J = l ; 

p—1 n k 

fc=l ii,...,ife+i=l (äi,...,(5fc)6AAP-i(fe) J = l 

p—1 n k k 

E E D^f^ E E^^^=+M^'^= n^s'^^. 

fc=l ii,...,ifc = l ((5i,...,<5fc)eA^P-i(fc) s=l i=l, 

/3:=eijH he;^ JT^s 
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For A; = 1, . . . ,p we set 

n k 



n:= E n^%-^^. (10) 



Now Our goal is to prove that: 

j^D,Vi = pTk (11) 

Our strategy of proof is as follows. We will define Si,. . . ,Sp, such that 

d ^ 

-Da^i = Y,Sk, S,=T„ i = l,...,p. (12) 
fe=i 

We set 

n k k 

^i = E E D^f^ E E^^-.^+MVi^n^s-^*.- 

fc=l ii,...,jfc=l (5i,...,<5fc)eAAp-i(fc) s=l j=l, 

ß~e-ii-\ l-e»j. JT^s 

n k 

Sp= J2 E ß'^/^-ß(c)^..+. E n^s^^r 

/s=p-l ii,...,ifc+i = l (5i,...,Äfc)eA/'P-i(fe) j=l 

For m = 2, 3, . . . , p — 1 we set 

n k 



fe=m-l ii,...,ifc+i = l (Äi,...,5fe)eArP-i(fe) j=l 

+ E E D^fi E E^f.^+w'^'. n^s-^^.- 

k=m ii,...,ifc=l (5i,...,5fc)eArp-i(fc) s=l j=l, 

ß-—ei^-\ hei^, j^s 

It remains to show that Si = Ti iov i = 1, . . . , p. Consider first i = 1. Recall 
that jVP-i(l) = {(l,2,...,j9- 1)}, hence 



s=l s=l 

Therefore 

Si=Ti. (13) 
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Consider now i = p. For an arbitrary s > A/'*(s) contains only one element 
((1), (2), . . . , (s)). Thercfore we obtain 

n p—1 

n, ■■■,«?=! (äi,...,(5p_i)e7Vp-i(p-i) j=i 

n p—1 
ii,...,ip = l j = l 

Since a = 6 + (c), where c = (op), hence 

n p 

ii,...,ip = l j = l 

n p 

= E ^"^+-+"''/.: E I[Das,^^,-T, 

il.---.ip = l (5i,...,äp)eAAP(p) j=l 

Consider now m = 2,3, . . . ,p — 1. We have 

n m— 1 

= l (5i,...,(5„_i)e7VP-i(m-l) J = l 

n m m 

+ E ^'^-+-+'^'"/. E E^^^,+w'^»= n^^v'., 

= l ((5i,...,(5„)eAAp-i(m) s=l J = l, 

Using decomposition J\fP{m) = A U i? as in ([5]) we obtain 

n m 

E ^"^^•••■'""/. E n^%-^^. 

n m 

+ E ^^-+-+^""/» E n^-.^». 

n m 
ii,...,i„, = l (i5i,...,(5„,)eAAP(m) J = l 

We have shown that = 5"^ for z = 1, . . . ,p. This finishes the proof. i 



4 C^-Lohner algorithm 

4.1 Why one needs an C^-algorithm? 

There are several effective algorithms for the computation of rigorous bounds 
for Solutions of ordinary differential equations, including Lohncr nicthod [Loj . 
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Hermite-Obreschkoff algorithm |NJ] or Taylor models [BM] , For C-computa- 

tions the number of equations to solve is equal to n I 1 hence, even for 

r = 1 direct application of such an algorithms to equations for variations (|14p 
leads to Integration in high dimensional space and is usuaUy inefficient. Let us 
recall after [Zl See. 6] the basic reason for this. In order to have a good control 
Over the expansion rate of the set of initial conditions during a time step these 
algorithms, while being C°, are 'internally'(or higher for Taylor models), 
because they solve non-rigorously equations for (|^) - the variational matrix of 
the flow. This effectively Squares the dimension of phase space of the equation 
and impacts heavily the computation time. But as it was observed in [Z] the 
equations for partial derivatives of the flow can be seen as non-autonomous and 
nonhomogenous linear equations, therefore we do not need additional equations 
for variations for them. As a result the dimension of the effective phase space 

/ ^ ^^ ^ \ 

for Our C'-algorithm is given by n ( ^ 1 and not a Square of this number. 

Another important aspect of the proposed algorithm is the faet that the 
Lohner-type control of the wrapping effeet is done separately for x-variables 
and variables Daip. This feature is not present in the blind application of C" 
algorithm to the System of variational equations and it turns out that this offen 
practically Switches off the control of the wrapping effeet on a;-variables, as 
various choices used in this control become dominated by the variables. 

In [Z] a C^-algorithm has been proposed. Here we present an algorithm for 
computation of higher order partial derivatives. 

4.2 An outline of the algorithm 

Let US fix r < K and considcr the foUowing System of differential equations 

d n k (14) 

k=l iu...,ik=l (<5i,...,<5fc)eA"<'(fc) J = l 

for all a e J\f^\ d = 1, . . . , r. 

Our goal is to present an algorithm for Computing a rigorous bound for the 
Solution of with a set of initial conditions 

lDipiO,xo) -Id (15) 
[Daip{0,xo) =0, forae7V2"U...U7V;'. 

In the sequel we will use the foUowing notations: 

• if a Solution of System (fT4)) is defined for t > and somc xq E R", thcn 
for a e A/" by Va{t,XQ) we denote Da(p{t,xo) 
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• for [xq] C M" by [Va{t, [xq])] we will denote a set for which we have 
Va{t,[xo]) C [so])]. This set is obtained using an rigorous nume- 
rical routine described below. 

The C-Lohner algoritlim is a modification of C^-Lohner algorithm [Z . One 
Step of C-Lohner is a shift along the trajectory of the System (fH|) with the 
foUowing input and Output data 
Input data: 

• tfc - a current time, 

• hk - a. time step, 

• [xk\ C R", such that Lp{tu, [xq]) C [xk], 

• [Vk,a] = [Vk.a{tk,[^o])] C M", such that Daipitk,[xo]) C [Vk,a] for a e 
Afl'U ...UJV^. 

Output data: 

• tk+i = tk + hk - a new current time, 

• [xk+i] C M", such that (p{tk+i, [xq]) C [xk+i], 



[Vfc+l.a] = [Vk+l.aitk+1, [Xo])] C R", such that Da^itk+1, [Xo]) C [Vfc+l.a] 

for a e A/;" U . . . U 7V;\ 



We will often skip the arguments of Vk^a when they are obvious from the context. 

The values of [a;fe+i] and [V^+i^a], a G JV" are computed using one step 
C^-Lohner algorithm. After it is done, we perform the foUowing Operations to 
compute [Vfc+i,a] for a e 7V2" U . . . U A/"" 

1. Find a rough enclosure for Daip{[0, hk], [xk])- 

2. Compute [Vfc+i,a], this will also involve some rearrangement computations 

to reduce the wrapping effect for V \Mo\ ILoj . 

5 Computation of a rough enclosure for Da^p 

For a fixed multipointer a E A/J Equation can be written as foUows 

^Daifiit, X) = Ba{t, X) + A{t, x)Daip{t, x) (16) 

where 

d n k 

k=2 ii,...,4fc = l (Äi,...,<5fc)eA/'<*(fc) J = l ^ > 

A^Dfoip 

The procedure for Computing the rough enclosure is based on the notion of 
a logarithmic norm, which we give below. 



10 



Definition 7 \HNWf For a Square matrix A the logarithmic norm ß{A) is de- 
fined as a limit 

\\Id + Ah\\-l 
H(A) = limsup 

where \\ ■ \\ is a given matrix norm. 

The formulas for the logarithmic norm of a real matrix in the most frequently 
used norms are (see |HNWj ) 

1. for ||a;||i = J2i ^^iA) = maxj{ajj + J2i^j \aij\) 



2. for m ||a;||2 = l^^ipj m(^) is equal to the largest eigenvalue of {A + 
A^)/2 

3. for ||a;||oo = max^ \xi\, ß{A) = maxi(aij + 

In Order to find bounds for Da^ we use the foUowing theorem [HNW) Thm. 
1.10.6] 

Theorem 8 Let x{t) he a Solution of a differential equation 

x'{t) ^ f{t,x{t)), xeM" (18) 
Let 1^(1) be a piecewise differentiable function with values in M". Assume that 

M V)] < m for V e [x{t),iy{t)] 



dx 



W{t)-f{t,vm<m, 



where by n{A), we denote a logarithmic norm of a Square matrix A G M"^". 
Then for t > to we haue 

\xit) - iy{t)\ < e^(*) (^\xito) - i^ito)\ + e-^(^)<5(s)ds^ , (19) 

with L{t) = J*^ l{T)dT. 

We apply the above theorem to Equation (fT6|) to obtain 



Lemma 9 Let us fix x £ M". Assume that |i3a(^, a;)| < 5{t) and ^{A{t,x)) < 
l{t), then for t > to 



\Da^{t,x)\ < |Da^(to,x)|e^(*) +e^(*) / e~'^^^^6{T)dT (20) 
with L{t) = //^ l{T)dT. 
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Proof: Consider Equation (|16p and a homogenous problem for 

^ f(t,w) -.^ A(t,x) -w, weM". (21) 

dt 

Using Theorem |S] we can estimate the difference between any Solution of (PT|) . 
w, and a Solution of denoted by ßa(/J. 

- wit)\ < \Da^{to) ~ z«(to)|e^(*) + e^(*) T e-^(-)<5(r)dr. (22) 

Jto 

After a Substitution w(t) — 0, which is a Solution of the homogenous equation, 
we obtain our assertion. I 

Usually, we do not have any control over the time dependence of S and l, 
hence we will use the following 

Lemma 10 Assume that \Ba{t,x)\ < S and ß{A{t,x)) < l for t e [0, /i] then 
for t Cz [0,h] we have 

\Daipit,x)\<\Daip{0,x)\ma^il,e''') + 6^-^, if l ^ 0, (23) 

or 

\Daipit,x)\ <\Da(piO,x)\+St, whenl = 0. (24) 



5.1 The procedure for the computation of the rough en- 
closure for V. 

The procedure for the Computing of the rough enclosure is iterative, which 
means that given a rough enclosure for (^([0, /ifc], [x^]) and rough enclosures 
Daip{[0,hk],[xk]) for all a G A/"" U . . . U Afp we are able to compute the rough 
enclosure for DafdO, hk], [xk]) for a € Afp+i- 

The procedures for computation of the rough enclosures of ip{[0,hk], [xk]) 
and Da95([0, /ife], [xk]) for a e A/"" has been given in [Zj. Below we present an 
algorithm for Computing [Ea] for a £ M2 U . . . U A/"" . 
Input Parameters: 

• hk - 8i time step, 

• [xk] C M" - the current value of x = ^p{tk, [xo]), 

• [Eq] C M" - a compact and convex such that (p{[Q, hk], [xk]) C [Eq] 

• [Ea] C R", a e TVi" U ... UTVp" such that DafHO, hk],[xk]) C [Ea] for 
aeAf^U...UAf;\ 

Output: 

• [Ea] C R", a e 7V;Vi such that 

Da^{[0,hk],[xk]) C [Ea] 
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Before we present an algorithm let us observe that for a fixed a € J^p+i, Ba 
defined in p?)) could be seen as a multivariate function of t, x and Vh = Di,(p for 
6 G A/fU. . More precisely, put rUp := jj {A/"" U . . . U AT^ }, where ji Stands 

for number of elements of a set. Recall that, we have defined by ([ü]) a linear 
Order in A/"". Hence, there is a unique sequence of multipointers bi,. .. ,bmp, 
such that bi e A/"" U . . . U Afp ior i — 1, . . . , m^, 61 < 62 < ■ • • < bmp and 6^ ^ &j 
for i ^ j. 

Let US define 



rrip + l 



by 



p+1 n 
fe=2 ii,...,ifc = l 

E n (-^, 

(5i,...,Äfe)6A^J'+i(fc) i=l 



(25) 



and 



Fa{t,X,Vb^,. . . ,U6„J = Bait,X,Vbi,. ..,Vb^)+ D f {ipit, x))Vait, x) (26) 

Algorithm: 

To compute [i^a] for a £ ■N'p+i we proceed as follows 

1. Find l > (max^gfßQ] ^i{Df{x))). 

2. Compute 5a > max||_ßa||, i-e. 

Öa > 



max 

{x,vtj^,...,vb^ )e[Eo]x[Eb^]x---x[Eb„ ] 



BaiO,X,Vb^,...,Vb,^ ) 



For example, if a = (j. c) e A/"^, then (5a should be such that 



öa > max 

2;e[Eo],-"i6[-E(i)],---,i'r.6[ß(„)] 



E 

r,s — 1 



9V 
dxrdx. 



3. Define [Ea]i — [— 1, l](5a^— for i = l,...,n, where [E'qJ^ denotes i-th 
coordinate of \Ea\- 

One can refine the obtained enclosure by 

{Ea\ ■■= {[0,hk]Fa (O, [Eo], [Eb,], . . . , [Eb^p])) n [Ea] 
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Indeed, for z = 1, . . . , 71, i G [0, hk] and xo G [Eq] we have 

Da(pi{t,xo) = ßa-^Ki^a^o) - -Da'^(0,a;o) 

= *(-F'a)i (6'i,a;o,^bi<^(öi,a;o), . . . , Db^^(p{9i,xo)) 

= *(-P'a)i {0,(p{Oi,XQ),Db^(p{di,xo), . . . ,Db^^(p{9i,xo)) 

for some G [0, C [0, hk]- In the above we have used the fact that 

Fa(t,X,Vi, . . . , Wmp) = Fa{0, ip{t,x),Vi, . . .^Vm^). 

Since ip{9i,Xo) G [£^0] and Di,^(p{0i, xq) G [£^f,J for j = 1, . . . , nip we get 

Da^^it,Xo) G [0,/lfc] (£^a), (o, [Eo], [Eb,], [Eb„j) 

6 Computation of [V^+i] 

6.1 Composition formulas 

For any p-times continuously differentiable functions /, g : R" R" and a G Afp 
we have 

p n k 

Da{fog) = J2 E E n^%5^. (27) 

fe=i ii,...,ifc=i (5i,...,äfe)e7Vp(fc) i=i 

We can apply the above forniula to / = (p{hk, ■) and g = 'y3(ifc, ■) to obtain 

p n 

Va{tk+hk,Xo) ^ VX-i(eij+...+eij(/lfc,a:fc) 

fe=l ii,...,ifc=l 

k 

(5i,...,Äfc)eA^p(fc) j=i 

for aU Xq G [xq]. Using notations [Vfe+i,a] := [K(ifc + /ifc, [a;o])] and [Vk,a] = 
[Va(tk, [xq])] we can rewrite the above equation as 

p n k 

[Vk+l,a]^Yl E ^/^-He,,+...+e,j{hkdxk]) II , . 

fc=i n,...,ifc=i ((5i....,Äfc)eA"p(fc) j=i 

(28) 

where A is defined by 

6.2 The procedure for computation of [V^+i] 

We introduce new parameters Od - the order of the Tayfor mcthod used in 
computations of Va for a G JV^ ■ It makes sense to take oi > 02 > • • • > . 
Input Parameters: 
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• /ifc - a timc Step, 

• [xk] C M" - the current value o{ x = ^{tk, [xq]), 

• [Vk,a] C K" - a current value of Vk,a{tk, [xq]), for a G U . . . U TV;" 

• [£'o] C M" compact and convex, such that >f{[0, hk], [xk]) C [Eq] - a rough 
enclosure for [xk], 

• [Ea] C K", compact and convex, such that Da(p{[0,hk], [xk]) C [Ea], for 
ae7V{^U...U7V;". 

Output: [Vfc+i.a] C M", such that 

Vaitk+hk,Xo) e[Vk+i,a] (29) 

for Xq e [xo] and a e TVi" U . . . U Af^ . 
Algorithm: We compute [Vfe+i] as foUows 

1. Computation of Va{hk, [xk]) using Tayfor method for Equatfon i.e. for 
a G Afp we compute 



i— 1 

+ (^l^^"(0,[i^o],[i^.J,...,[i^.,,_J). 

where Vt. = for bi e U . . . U TV^^i and Vq) = for j = 1, . . . , n. 
Observe that 

Vaihk,[xk])(l[Fa] (31) 

Indeed, using Tayfor series expansion we obtain that for Xk G [xk] and 
j — 1, . . . ,n holds 

j 2 1 

(K ), (/ife , xfe ) = ^ — ^ (F„ ), (0, xfe , Vfc, (0, xfe ) , . . . , _ ^ (0 , ) ) 

+-^^^^^(^a)j(ö„a:fc,14,(Ö„a:fc), . . . ,14„^_^(ö„Xfc)) 
for some e [0, hk]- Observe, that 

|^(F„),(ö., xk,n, {e..,xk), . . . , (ö., xfc)) 

= (0' ^(^^^' ^fc)' (O^^^k), . . . , 0, (ö„ xfc)) 

Using ip{9.i,Xk) G [i?o] and Vfc^(6'i,a:fe) G for s = l,...,mp_i we 

obtain our assertion. 



-1 
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2. The composition. Put 



[Jfe] :=([F(i)],...,[i^(„)]f 
Using (HB]) for a e J\f^ we have 

[Vk+l,a] = [aa] + [Jk] ■ [Vk,a] 

where 



(32) 



^J=E E [^A-(e,,+...+e.j] E nK^J,. (33) 



(äi,...,5fc)eAAp(fc) i=i 



In Our implementation of the algorithm we use the symbolic differentiation to 
obtain formulae for Daf. Ncxt, using the automatic differentiation we compute 
§iFa{t,x, Vb^{t,x), Vb^^_^ {t,x))\t=Q which appear in 



6.3 Rearrangement for Va - the evaluation of Equation 

It is weh know that a direct evaluation of Equation ([32]) leads to wrapping effect 
|Mo[ ILo] . To avoid it foUowing the work of Lohner [ Loj we will use the same 
scheme as it was proposed in [Zj. 

Namely, observe that Equation ([32]) has exactly the same strueture as the 
propagation equations for C^-method (see |Zl Section 3]). Moreover, all vectors 
Vk,a, for a G A/"" U . . . A/"" 'propagate' by the same [Jk] as did the variational 
part in [Z], hence it makes sense the same approach. 

To be more precise, each set [Vfc^a], for o G A/"" U . . . U A/"" is represented in 
the foUowing form 

[14 ,a\ — '^k.a + [Bk][rk,a] + Ck[qk,a] 

where [Bk] is interval matrix, Ck is point matrix, Vk,a is a point vector and 
'f'k,a,<lk,a are interval vectors. Observe that [B^] and Ck are independent of a. 
In the sequel we will drop index a. Equation (|32|) leads to 

[Vk+i] = H + [Jk\{vk + [Bk][rk] + Ck[qk]) (34) 

Let m([z]) denotes a center of an interval object, i.e. [z\ is interval vector or 
interval matrix and A([2;]) = [z] — m{[z\). 

Let [Q] be an interval matrix which contains an orthogonal matrix. Usually, 
[Q] is computed by the orthonormalisation of the columns of iin{[Jk])[Bk]. 

Let 

[Z] = m{[Jk])Ck 
Ck+i = m{[Z]) 
[Bk+i] = [Q] 
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Then we rearrange formula ([M]) as foUows 

(35) 



[s] = [a] + [Jk]vk + A{[Jk])[Vk] 
Vk+i = m{[s]) 

[n+l] = miA{[s]) + Ai[Zmq,]) + {[Q^]m{[.mB,])[r,] 



Summarizing, we can use the foUowing data structure to represent ip{tk, [xq]) 
and Daifiitk, [xa]), for a e Wi" U . . . U TV;" 

type CnSet = record 

vo,ro,qo: IntervalVector; 
Co, Bo,C, B : IntervalMatrix; 
{va,ra,qa ■ IntervalVector}aeAAi"u...uAA;' 

end; 

The set ip{tk, [xq]) is represented as vq + Borg + CqQq, the partial derivatives 
Dafitk, [xq]) are represented as Va + Bra + Cqa- The matrices B, C are common 
for all partial derivatives. 

Notice, that if we start the C computation with an initial condition (|15p then 
there is no Lipschitz part at the beginning for the partial derivatives. Hence, 
the initial values for C and B are set to the identity matrix and the initial values 
for qa,ra are set to zero. 

If the interval vectors become 'thick' (i.e. theirs diameters are larger than 
some threshold value) we can set a new Lipschitz part in our representation (it 
must be done simultaneously for all Da(p) and reset in the foUowing way 

qa - ra + {B^C)qa, forae AA{^U...UAC 
ra = 0, forae A/'f 
C = B 
B = Id 

A similar change of the Lipshitz part may be done when vectors become thick 
in comparison to qa- 

7 Derivatives of Poincare map 

Consider a differential equation 

x'^fix), xeR", /eC^''+^ (36) 

Let : M X M" ^ R" be a (local) dynamical System induced by Let 
a : K" ^ R be C^-map. Put U = {x \ a{x) = C}. 

Definition 11 We will say that H is a local section for the vector field f at 
yoeu if 

(Va(2/o)|/(2/o)) ^0. (37) 
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Assume xq € and S K are such that II is a local section at ip{to, xq). 
Considcr an implicit equation 

a{ip{tp{x),x))=C. (38) 

It foUows easily from ([57]) and from the implicit function theorem that there 
exists a uniquely defined tp : R"-6->R in a neighborhood of xq, such that 
tpixo) — to- The function tp is as smooth as the flow ip. We will refer to 
tp as to the Poincare return time to section II. 

We define a Poincare map P : R" D dom (tp) — » R" by 

P{x)=^{tp{x),x). (39) 

Usually the Poincare map is defined as a map P : Iii -«^112, where 111,112 are 
local sections in R". The approach taken here, i.e. treating the Poincare map 
as map P : R"-e-^R" allows us to not to worry about the coordinates on local 
section. 

In this section we are interested in the partial derivatives of P defined by 
(EU). 

From (|39p wo can compute ||i and we obtain 

^{x) = f,iPix))^ix) + ^{tp{x),x). (40) 



We need We differentiate §SI to obtain 




7.1 Higher order derivatives of the Poincare map 

To make formulas transparent we will drop arguments of functions in this sec- 
tion, but reader should be aware that for tp and its partial derivatives the 
argument is x, for and Daf the argument is always the pair {tp{x),x). 
From (PÜ)) we obtain 

d d 
Dü,c)P ^ Q^'PD{3)tpD(^^^tp + —D(^^)ipD^j)tp + —ipD(^j,^-^tp 

d 
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It is easy to see that partial derivatives of high order give rise to quite complex 
expressions and it is not entirely obvious how to organize it in some coherent 
and programmable way. For this purpose we use the following 

Lemma 12 For a multipointer a G Mp we have 

DaP = Daip+pDatp 

+ Sfe=2 "ät^ E((5i,...,(5fc)GAAp(fc) nj=l -C'oä^ (43) 
+ J2k=2 J2{Si,...,Sk)eN'Pik) J2s=l dt^-i^as^y^Hj^s ^as.tp 

Proof: By induction onp. Forp = 1 formula (|43p is equivalent to ([^0]) . because 
thc two last sums are taken over empty set. Assume p3)) holds true for some 
p>l and fix a G J^p+i- Om goal is to show that 



DaP = i?i + i?2 + i?3 



where 



d 

Rl = Da(p + -;-ipDatp 

ot 

p+1 ^fc fc 

k=2 (äi,...,5fc)eAAp+i(fc)i=i 
P+i k Qk-i 

^3 = E E Eä^^-.^n^-.^^ 

Write a — ß + j, where ß G Afp and 7 — (op+i) G A/"". From the induction 
assumption we have 

DaP - D^{Dß^+§^Dßtp) 

+ ^1 {YTk=2 ^'^E(5i,....äfc)eAAp(fc) 11^=1 Dß,.tp^ 
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where 



^1 

^3 
S7 



§-^DßCpD^tp 



ö^^ipDßtpDytp 
■ß^D^ipDßtp 

TJk=2 m^^DjtpJ2iSi,....Sk)eN-p{k) 11^=1 Dß,.tp 



Efc=2 m^'fiT,{Si.,...,Sk)eMp(k) T,s=i Dßs,+fipl\j=i Dß.tp 




Obviously Ri ^ Si. We will show that R2 = Sj, + Sq + Sj and R3 = S2 + S4 + 



Denote by i = 2, 3 a part of sum Ri with fixed k — 2, . . . ,p + 1. 

Similarly, let us denote by Si^k a part of sum 5*^, i = 5, . . . , 10, for k ~ 2, . . . ,p. 

Using decomposition of A/'^"^^(2) as in ^ we obtain that i?2,2 — S3 + 87,2- 
Similarly, using ([8]) we observe that i?2,fe = <5'6,fc-i + Sr,k for k = 3,...,p. 
Finally, since AfP+^{p + 1) = {((1), (2), ...,{p+ 1))} and 7 = (op+i) we find 
that i?2,p+i = Sß.p. This shows that R2 = S3 + Se + Sj. 

It remains to show that R3 ^ S2 + S4 + S5 + Ss + Sg + Sw . We will classify 
possible terms by the fact, where p + 1 appears in Si, i — 1, . . . ,k and how this 
Si enters in R3 as Ss or Sj . There are four cases 

1. 6s = + 

2. S, = {p+l) 

3. p+leS„ \Ss\>2 

4. p+le Sj, \öj\ > 2 

Let US fix /c = 2. Let {61,52) G AfP'^^{2). The term for case 1 is 5*4, for case 2 
is S2, case 3 is ^8.2 and case 4 is S'io,2- Hence, i?3^2 = £"2 + <S'4 + 5*8,2 + 5'io,2- 

For k — 3, . . . ,p and fixed {61, ... ,6k) & N'P~^^{k) we have: case 1 is given 
by S's.fc-i, case 2 by Sg^k-i, case 3 by Ss.k and case 4 by 5io,fc Hence, for 
fc = 3, . . . ,p we have i?3,fc = 55,fc_i + Sg^k-i + Ss^k + S'io./c- 

Finally, for fc = p + 1 we observe, that i?3.p+i = S^^p + Sg^p. Indeed, in this 
case {61, . . . , 6p+i) — ((1), (2), . . . ,{p + 1)). Hence, either for (5^ = 7 we have 
term S^^p and 6s j we have Sg^p. 

We have showed that R3 = 5*2 + ^4 + ^5 + ^8 + 6*9 + 5*10 and the proof is 
finished. I 



S5 + Ss + Sg + Sio 
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Hence, if we know all the partial derivatives of tp up Order p we can compute 
the partial derivatives of the Poincare map up the same Order. In next subsection 
we show how to compute partial derivatives of ip for affine sections. 



7.2 Partial derivatives of for affine sections 

Assume a : R" — s- R is an affine map given by 

n 

a{x) = ao + ctiXi. 
1=1 

This is a quite restrictive assumption about sections, but it leads to relatively 
simple formulas for Datp and it is sufficicnt for the applications wc havc in 
mind. 

Lemma 13 For a multipointer a G Mp holds 

+ EL2 (V"l^<^) E(Äi,...,ä.)eA^P(fc) 11^ = 1 Da.^tp 
E'=l (Va||^öa,.^)n,^.öa./P 

Proof: The proof is a direct consequence of Lemma [T^] and ([55]) . Since a is 
afhne, by differentiating of a{P{x)) = C we get (ValDaP) — 0. Using formula 
(|43l) for DaP we obtain our assertion. I 



Fix [x] C M" and assume we have a rigorous bound for tp{[x]) E [ti, 12] (see 
[Z| Section 6] for more details on this). Lemmas fT3l and fT2l show that given rigo- 
rous bounds for the partial derivatives Daip{[ti, t2], [x]) and -ß^Daipdti, t2], [2^]) 
up to some order p we can compute recursively rigorous bounds for the partial 
derivatives of tp{[x]) and up to the same order. Notice, that -^Da^p are 

given by Taylor coefficients of the Solution of (|14p with initial conditions P([a;]) 
for C*^ part and Da^{tp{x), [x\) for equations for variations. Hence, thcse coef- 
ficients can be easily computed using the automatic differentiation algorithm. 



8 Applications. 

One of the typical invariant sets in hamiltonian mcchanics are invariant fori. 
However, the existence of invariant torus in a given System is offen difficult to 
prove despite the fact that the theory is quite well developed. Probably the 
best work in this dircction was done by Celletti and Chercia [CCH ICC2| , where 
the an effective application (computer assisted proof) of KAM theory to the 
restricted three body problem modelling System consisting of Sun, Jupiter and 
asteroid 12 Victoria was given. Our aim here is more modest as we locus on the 
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invariant tori emanating from the elliptic fixed point satisfying siiitable twist 
condition. 

In this section we show that the rigorous computations of partial derivatives 
of a dynamical System up to Order 3 or 5 can be used to prove that in a particular 
System an invariant torus exists around some elhptic periodic orbits. In this 
section this will be done for the forced pendulum equation and the Michelson 
System. 

8.1 Area preserving maps on the plane, normal forms and 
KAM theorem 

Definition 14 Let f : M.^ ^ M.^ be a smooth area preserving map, such that 
f{p) — P- Let X and fi he eigenvalues ofdf{p). Following JSM/ we will call the 
point p 

• hyperbolic i/ A, /i G M and X ^ fi, 

• elliptic if X =JI and A 7^ /i, 

• parabolic if X — fi. 

The following KAM theorem will be the main tool to prove the existence of 
invariant tori in this paper. 

Theorem 15 \SM\ %32] Consider an analytic area preserving map / : M'^ — > 
R^, f(r,s) = (ri,si) where 

ri = rcosa — ssina + 02;+2 

Sl — r sin a + s cos q; + 02;+2 (44) 
/ 

E/ 2 
Jk[r + s ) 

k=0 

and O21+2 denotes convergent power series in r, s with terms of order greater 
than 21 + 1, only. 

If at least one 0/ 71 , . . . , 7/ is not zero then the origin is a stable fixed point 
for map f . Moreover, in any neighborhood U of point there exists an invariant 
curve for map f around the origin contained in U . 

The next theorem and its proof teils how to bring a planar area preserving 
map in the neighborhood of an elliptic fixed point into the form (j44p . 

Theorem 16 JSM[ %23] Consider an analytic area preserving map / : — )■ 
such that /(O) = 0. Let A, A he complex eigenvalues of Df{0), such that \X\ = 
|Ä| = 1. If X^ ^ \ for k — 1, . . . , 2Z + 2, then there is an analytic area preserving 
Substitution such that in the new coordinates mapping f has form (|44p . 

The proof of the above theorem is constructive, i.e. given the power series 
for / at an elliptic fixed point one can construct explicitly an area preserving 
Substitution and compute the coefScients 70, . . . , 7; in ([44]) . An explicit formula 
for the coefficient 71 in the above normal form is given in Appendix [X] 
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8.2 The existence of invariant tori in forced pendulum. 

Consider an equation 

e^- sm{9) + sin(cji) (45) 

Observe that pS)) is hamiltonian. 

Let US denote by : R^-e->R^ the Poincare map for Equation (|45| with 
a Parameter uj, i.e. — (p{2Tr/LLj, •), where (/j : R x R^-e-^R^ is a local flow 
induced by flS]) . Observe that (|l5)) is nonautonomous, but it is equivalent to 
first Order System of autonomous ODE given by 

V 

- sm{e) + sin(wt) (46) 
1. 

In the sequel all rigorous computations for ((45|) will be in fact performed for the 
System ([iS]). 

Observe that to any invariant closcd curve for corresponds and invariant 
2-torus for (|45)) . 

Consider a set of parameter values 

rii = [2,2.994], = [3,3.997], r!3 = [4,8] 

The foUowing lemma was proved with Computer assistance 

Lemma 17 For all parameter values w g ri there exists an elliptic fixed point 
Xuj G R^ for P^. Moreover, there exists an area-preserving Substitution such 
that in the new coordinates the map fuj{x) — Puj{x + x^j) — Xuj has the form (|44p 
with l — 1 and 71 7^ 0. 

Before we give the proof, let us briefly comment about the choice of the parame- 
ter set $7. For parameter values slightly lower than 2 we observe the parabolic 
case, i.e. there exists a parameter value uoi for which eigenvalues of the deriva- 
tive of are equal to — 1. In two gaps in below 3 and 4 we have resonances 
of low Order. Namely, we have parameter values with an elliptic fixed with ei- 
genvalues to e^^^/^ = ^ ± ^i and e*"/^ = ±i, respectively. Clearly, in a 
Computer assisted proof we need to exclude a small interval around those Para- 
meters. For > 4 it seems that the interval f^a can be extended much further 
to the right without any difflculty 

Proof of Lemma ll7t A Computer assisted proof consists of the foUowing steps. 
We Cover the set 51 by 9910 nonequal subintervals oji. Diameters of üj^'s were 
relatively large for values far away from the parabolic cases and very small close 
to them. For a fixed subinterval LOi we proceed as follows 



de 

ds 
dv 

ds 
dt 

ds 
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1. Let üj denote an approximate center of the interval Wi. We find an approxi- 
mate fixed point for P^d using the Standard nonrigorous Newton method. 
Let US denote such a point by Xi. 

2. We define a box centered at Xi, i.e we set Vi := Xi + [— £^,6^]^, where 

> depends on subinterval uji - the values we used are from the interval 
[5 • 10~^,3 • 10"'^], depending on whether Xi close to parameter values 
corresponding to parabolic cases. 

3. Using the C^-Lohner algorithm we compute the Interval Newton Operator 
[MöllNlE] Ni := N{P^^ -ld,Xi,Vi) and verify that N, C intw^. This 
proves that for all lo G uji there exists a unique fixed point e Ni for 

4. Using the C'^-Lohner algorithm we compute a rigorous bound for P^. (Ni) 
and D'^P^^ (iV^), a G UNj UN§. Hence, we obtain a rigorous bound for 
the coefficients in 



5. We show that an arbitrary matrix M E DP^.{Ni) has a pair of complex 
eigenvalues A, A which satisfy A'^ ^ 1 for fc = 1, . . . ,4. From Theorem [Tül 
it foUows there exists an area-preserving Substitution such that in the new 
coordinates the map for u; G uii has the form ([44]) with ^ = 1. 

6. We compute a rigorous bound for 70 and 71 which appear in the formula 

and verify that for co E oji holds 71 7^ 0. 

The rigorous bounds for the values of 71 on are 



A Computer assisted proof of the above took approximately 95 minutes on the 



As a straightforward consequence of Lemma [17] and Theorem [15] we obtain 

Theorem 18 For all parameter values lü £ there exists an elliptic fixed point 
x^ G for P^ . Moreover, any neighhorhood of point x^ contains an invariant 
curve for P^^ around x^j ■ 




7i(r2i) c [0.29930416771330087,30.118260918229566] 
71 (172) C [0.099747909112924596,0.56550301088840627] 
71(173) c [0.18574835001593507,0.4129279974577012] 



Pentium IV 3GHz processor. 
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8.3 Higher order normal forms. 

In the previous section it was shown that computations are sufficient to 
prove that for Ij45p a family of invariant tori exists. However, it may happen 
that the coeiBcient 71 in the normal form vanishes. In this Situation we may try 
to compute higher order normal form. As an cxamplc wo consider a pcndulum 
with a different forcing term, 



Theorem 19 Let he the Poincare map for (j47p . For all parameter values 
w € il* = [2.9957694795, 2.9957694796] there exists an elliptic fixed point Xuj G 
for P^ . Moreover, any neighbourhood of point Xuj contains an invariant 
curve for P^ around x^ . 

Proof: The main concept of the proof is the same as in Lemma [TTI Using the 
nonrigorous Newton method we find an approximate fixed point 

X = (-7.7491573604896152- 10^1^-0.54723831527031352). 

We set V = X + 3 ■ 10^^([— 1, 1] x [—1, 1]). Using the C^-Lohner algorithm we 
compute the Interval Newton Operator of P,^ — Id on w and we obtain that for 
allco eü^, N = N{P^ - Id, center(w), v) C (iVi, iVa), where 

Ni = [-5.1582932672798325,5.1582631625020222] • 10^1" 
N2 = [-0.54723831580217108,-0.54723831470891193] 

Since N d v we conclude that for all w G fJ* there exists a unique fixed point 
Xlj E N for the Poincare map. 

Using C^-Lohner algorithm we compute a rigorous bound for Pq^ (N) and 
D°'Pn^{N), a € U . . . U N5. Hence, we obtain a rigorous bound for the 
coefficients in 



We show that an arbitrary matrix M £ DPn, (N) has a pair of complex ei- 
genvalues A, Ä which satisfy A*"' 7^ 1 for fc = 1, . . . , 6. From Theorem [Tülit foUows 
there exists an area-preserving Substitution such that in the new coordinates the 
map fi^ for w G $7* has the form (|44| with l = 2. 

Next, we compute a rigorous bound for 71 and 72 which appear in the formula 
(|44l) and we get 



e = - sin{9) + sin{Lüt) + sin(2wi). 



(47) 




71 (O*) C [-5.3924276719042241,5.381714805052106] • 10 
72(f7*) C [199.95180660157078,199.99104965939162] 



1-6 



Since for uj E ft^,, 72 (ti») 7^ the assertion foUows from Theorem [T5l 
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The main Observation which makes this examplc interesting is that there 
exists (jj* S 51* for which 71 (w*) = and we cannot conclude the existence of 
invariant tori for all G 51* from computations. To be more precise, we 
computed the coefficient 71 for the parameter values ui = min fi* and lü2 = 
maxfi* and we get 

7i(t^i) € [-2.3559594437885885,-1.3593457220363871] • 10"** 
71(^2) e [2.9671154858524365 • 10"^ 1.2819312939263052 • 10"^] 

Since 71 exists for all G fl* and depends continuously on uj we conclude, that 
71 (w*) = for some lj* G 17*. 

8.4 Application to the Michelson System 

The existence of an invariant curve for a planar map / : ^ can be proven 
without assumption that / is measure preserving. The key assumption in the 
proof given in [SM is that any curve 7 around an elliptic point intersect its 
Image under /, i.e. /(7) H 7 7^ 0. Such a Situation is also observed in reversible 
planar map around an Symmetrie elliptic fixed points. 

Definition 20 An invertible transformation M : — > Q is called a reversing 
symmetry of a local dynamical System (j) : T x Q — > Cl. T = M or T = Z, if the 
following conditions are satisfied 

1. if {t, x) G dom (0) then {~t, S{x)) G dom (</)). 

2. S{<l>{t,x)) = <l>{-t,S{x))) 

Remark 21 In the discrete time case, the above two conditions are equivalent 
to identity 

M o f = f-^ o M. 
where f — (f)(l, •) is a generator of cß. 

Definition 22 Lei cj) : T xQ Q be a local (discrete or continuous) dynamical 
System. For a; G 51 put 

I{x) = {i G T : (i, x) G dom ((/))} 
0{x) = {(l){t,x) (^ü:t I{x)} 

The set O {x) will be called a trajectory of a point x. 

Definition 23 Assume S is an reversing symmetry for : T x O — > 51. An 
orbit O (x) is called S -Symmetrie orbit if O (x) — S{x). 

Remark 24 [La] In continuous case the orbit O (x) is S -Symmetrie if it contains 
a point from the set Fix(S') — {y : S{y) — y}. 
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Remark 25 JWi^ Lern. 3. 3] It is easy to see that ifQcflisa Poincare section 
for a R-reversible flow : R x — > f2 such that Q — R{Q) then the Poincare 
map P : O — > G is R\q -reversible. 

As we observed at the beginning of this section, an i?-reversible planar map 
may admit an invariant curve around an i?-symnietric elliptic fixed point. In 
reversible case a planar map admits the same normal form around Symmetrie, el- 
liptic fixed point as in the area-preserving case and the Substitution which tends 
the map to the normal form is exactly the same as we described in Appendix [Xl 
-for details see [S^ lBHS] , 

Consider an ODE 

X = y 

V = z (48) 

9 1 ? 

z — c — y ~ 

On one hand, the System is an equation for the steady State Solution of 
one-dimensional Kuramoto-Sivashinsky PDE and it is known in the literature 
as the Michelson svstem ^Mi) . On the othcr hand, this System appears as a part 
of the limit family of the unfolding of the nilpotent singularity of codimcnsion 
three (see |DIKlj ). 

The System (|48p is reversible with respect to the symmetry 

R: {x,y,z,t) ^ {-x,y,-z,-t) (49) 

and since the divergence vanishes it is also volume preserving. 

A dynamical System induced by (|48p exhibits several types of dynamics for 
different values of parameter. For sufficiently large c there is a simple invariant 
set consisting of two equilibria (±c\/2, 0, 0) and heteroclinic orbit between them 
|MC| . Lau |Lau| numerically observed that when the parameter c decreases a 
cascade of cocoon bifurcations occurs and at the limit value c « 1.266232337 
a periodic orbit is born through a saddle-node bifurcation. This hypothesis 
has been proved in [KWZj . The Computer assisted proof of this fact given in 
[KWZ] uses the algorithm presented in this paper in order to compute partial 
derivatives up to second order for a certain Poincare map. 

For the parameter value equal to one and slightly smaller than one it was 
proven m |DIK2[ IWin IWi2l IWi3] that the System has rieh and complicated 
dynamics including symbolic dynamics, heteroclinic Solutions, Shilnikov homo- 
clinic Solutions. 

However, as the bifurcations diagram presented by Michelson suggests jMii 
Fig.l] for all parameter values c G (0,0.3195) there are at least two elliptic 
periodic orbits with large invariant Islands around them. In this section we 
present a proof that such Islands exist for some ränge of parameter values. The 
main idea of the proof is almost the same as in the previous section. There are 
two main differences. First, the Poincare map will not be a time shift. Therefore 
computations of the partial derivatives of the Poincare map require Lemma [T^] 
and Lemma 1131 Second difference is: we use the shooting method instead of 
the interval Newton method for the proof of the existence of Symmetrie periodic 
orbit. 
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The aim of this scction is to prove the foUowing 

Theorem 26 For all parameter values from the set 

C = Ci U C2 = [0.1, 0.225] U [0.226, 0.25] 

there exists a Symmetrie elliptic periodic orbit for the Michelson System (j48p . 
Moreover, each neighhourhood of such an orbit contains a 2D tori invariant 
under the flow generated by the Michelson System. 

Let US define the Poincare section II := {{0,z,y) : z,y £ R}. Let Pc = 
(Pi, P2) : n-6^n be the Poincare map for the System with the parameter value 
c. Notice, that Pc is in fact a half Poincare map, which means that the trajectory 
of X crosses II in opposite directions when passing through x and Pdx), and 
therefore periodic orbits for the Michelson System corresponds to periodic points 
for P2_ 

Since the section II is invariant under symmetry {x,y,z) — s- {—x,y,—z), 
from Remark [21] the Poincare map is also reversible with respect to an involu- 
tion R{y, z) — {y, ~z). We will use the same letter R to denote the reversing 
symmetry of the Poincare map and the Michelson System. 

Let US comment about the choice of the set C. In the gap between intervals 
Ci and C2 there is a parameter value for which the eigenvalues of the Poincare 
map P^^ are ±z. Apparently at this parameter value we have a bifurcation and 
four periodic islands are born as it is shown in FiglT]- see also a movie mpp.mov 
available at [Wi4| which presents an animation of the phase portrait of Pc for 
the parameter values from the ränge [0.1, 0.25]. 

Proof of Theorem 1261 The main concept of the proof is quite similar to the 
one presented in Lemma[T71 We divide the set C of parameter values onto 20800 
nonequal parts (smaller when close to the bifurcation parameter c* and close to 
0.1 and 0.25). For a fixed subinterval Ci from the grid we proceed as foUows 

1 . Let c denote a center of the interval Ci . We find an approximate fixed point 
of P| using the Standard nonrigorous Newton method. Let us denote this 
point by {yi,Zi). 

2. Since the map Pc is reversible one can prove the existence of the fixed 
point for P^ using the shooting method as foUows. 

Let Fix(P) = {{y,z) £ U : R{y,z) = {y,z)} = {{y,0) e U : y e R}. 
Since Pc satisfies (Pc o P)^ = Id whenever the left side is defined, one 
can see that ii x £ Fix(P) and Pc{x) £ Fix(P) then P^ix) — x. Let us 
remark, that we always get an approximate fixed points {yi, Zi) resulting 
from the nonrigorous Newton method very close to Fix(P). We define 
two points ui = (jji — Ei, 0), U2 = {yi + £i, 0) £ Fix(P), where Si is a small 
number depending on ct and we show that Trz{Pci{ui)) ■ (Pc- (1x2)) < 0, 
where tTz is a projection onto z coordinate. Hence, if the P^ is defined 
on the set Ni — (0, [yi — ei,yi+ e«], 0) then for all parameter values c£ Ci 
there is a point Uc £ N which satisfies TTziPciuc)) = and therefore 
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Figure 1: Phase portrait of thc Poincare map Pc (top) before bifurcation for 
c = 0.225 and (bottom) after bifurcation for c — 0.226 with four periodic islands. 
Between those parameters resonant case occurs with eigenvalues equal to ±z. 
See also auxihary material |Wi4| . 
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Pc(uc) G Fix(i?). This shows that for all c G Ci there exists a fixed point 
for inside Ni provided Pc is defined on Ni, which will be discussed 
below. 

3. Using C'^-Lohner algorithm wc computc rigorous bounds for P^,{Ni) and 
D^'PKNi) for a e U U N^. This iraplics also that Ni C domPc.- 

4. We show that an arbitrary matrix M £ DPc-{Ni) has a pair of complex 
eigenvalues A, A which satisfy A'^ 7^ 1 for /c = 1, . . . ,4. From Theorem [Tül 
it foUows there exists an area-preserving Substitution such that in the new 
coordinates the map Pc for c e q has the form ((44|) with ^ = 1. 

5. We compute a rigorous bound for 70 and 71 which appear in the formula 
(|44| and verify that for c G q holds 71 ^ 0. 

The rigorous bounds for the values of 71 on C are 

7i(Ci) c [0.014515898754816965,157.76639522562903] 
7i(C2) C [1.1002393483255526,151.35147664498677] 

The Computer assisted proof of the above took approximately 7 hours and 50 
minutes on the Pentium IV 3GHz processor. 1 

9 Implementation notes. 

AU the algorithms presented in this paper have been implemented in C++ by 
authors and are part of the CAPD library [CAPD] . In particular, the package 
implements the computation of partial derivatives of a flow with respect to 
initial condition, partial derivatives of Poincare maps for linear sections and 
computations of normal forms for planar maps up to Order 5. 

The Implementation combines the automatic and symbolic differentiation in 
Order to generale a coefficients in Taylor series for the Solutions of the System 

m- 

Our tests shows that without difficulty we can compute partial derivatives 
up to Order 3 for an equation in 8-dimensional phase Space (which gives 1320 
equations to solve) on a Computer with 512MB memory. However, our current 
Implementation is optimized for lower dimensional problems. All the trees which 
represent formulas p4p are stored in the memory of a Computer. This speeds 
up computations because we do not need to recompute all the multiindices, 
multipointers and submultipointers in each step of the algorithm. Unfortunately, 
such an Implementation is memory-consuming. Therefore, higher dimensional 
Problems require a Computer with huge memory even for or computations. 
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A Explicit formulas for third order normal forms 
for a planar map 

The goal of this section is to give some details about the proof of Theorem [TÜl 
We want to present some formulas to give the reader the feehng about the 
necessary computations. 

Throughout this section we assume that the assumptions of Theorem llOl are 
satisfied. In the neighbourhood of / is given by a real, convergent power series 

f{x,y) = {xi,yi) 



Xi 



^^ük-uix'-y 



k=l 1=0 
oo k 

k=l 1=0 



Denote also by / : ^ a complex extension of /. Let A, Ä G C be 
complex eigenvalues of Df{0) and v,v € C corresponding eigenvectors (here 
bar denotes the complex conjugation). Then, using a linear Substitution of the 
form L = [v'^jv'^], we can change the coordinate System such that in the new 
coordinates the mapping / has the form 

f{^,v) = (Ae+p(C,?y),Ä77 + g(e,?7)) 



l^k-l 



k=2 1=0 
oo k 

k=2 1=0 

= qj,t fori,j>0. 

The last condition is a consequence of the invariance of C under the 
complex map /. We will refer to it as the reality condition. Namely, the set 
R2 c in the new coordinates ij) is given by ^ = 77 and the condition 
/(M^) C expressed in coordinates (^,77) is equivalent to ([5D|) . 

Assume now, that X'' ^ 1 for fc = 1, . . . , 4. Then an analytic area-preserving 
Substitution satisfying reality condition ([5D)) 

mz,v),-f{z,v)) = {zi,Vi) 

3 fc 

zi = z + ^^cßi^k-iz^v''-^ + ■■■ 

k=l 1=0 
3 k 

Vi = w + ^^V;,fc-iz't'''"' + ••• 

k=l 1=0 
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where 



-02,0 = (, 


^0,2 = 


-A%,2(A^ - 1)"' 


01,1 = (, 


ii,i = 


-pi,i(A-l)-i 


00,2 = (, 


!!'2,0 = 


P2,o(A2 - A)-i 


03,0 = (, 


^■0,3 = 


-A^ (po,3 + Pi,i0o,2 + 2go,20o,2) (A"* - 1)^^ 


02,1 = (, 


^1,2 = 


^2 i (Pl,2 + 2p2,O0O,2 + Pl, 101,1 + Pl, 100,2 


01,2 = <; 


^2,1 = 


-02,000,2 + 00,202,0 


00,3 = (, 


*3,0 = 


(P3,0 + 2p2,O02,O +Pl,102,o) (A^ - A)"^ 


brings / = 


= (/l,/2) 


to the normal form 



(z, v) — > {z{ao + a2zv), v{ßo + ß2zv)) + 0{{zvY) 

with 

ß2 ^ a2 = qi,2 + 2^2,000,2 + (71,101,1 + (71,100,2 + 2go,20i,i 
Finally, let 70 e M be such that A = ao = e*'''" and we compute coefficient 71 by 

-ia2 iß2 

71 = = -5- 

Qo Po 

From the proof given in [SM| it foUows that 71 <E K and the mapping / in 
coordinates (z,v) has the form 



f{z,v) = (^ze*('>'«+^i^''\i;e-'('^»+'^i^")) 



O4 



where O4 is a convergent power series with the terms of degree at least 4. 

Again, the coefRcients of f{z,v) satisfy reality condition ([5D|) . In order to 
express this normal form in terms of real variables we make a linear Substitution 

z ^ r + is, V — r — is 

and we obtain the normal form for / 

/(r, s) = (ri,si) + 04 

n = r cos(7o + 7i(r^ + s^)) — ssin(7o + 7i(r^ + s^)) 
Si = rsin(7o + 7i(r^ + s^)) + scos(7o + 7i(r^ + s^)) 



which agrees with 

The formulas for higher order terms 4>i,j 1 ^i.j (and for 72, which are not given 
here) has been computed in Mathematica. 
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